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Abstract 



After being considered as a nuisance to be filtered out, it became recently clear that biochemical noise 
plays a complex role, often fully functional, for a genetic network. The influence of intrinsic and extrinsic 
noises on genetic networks has intensively been investigated in last ten years, though contributions on 
the co-presence of both are sparse. Extrinsic noise is usually modeled as an unbounded white or colored 
gaussian stochastic process, even though realistic stochastic perturbations are clearly bounded. In this 
C) • paper we consider Gillespie-like stochastic models of nonlinear networks, i.e. the intrinsic noise, where the 

model jump rates are affected by colored bounded extrinsic noises synthesized by a suitable biochemical 
state-dependent Langevin system. These systems are described by a master equation, and a simulation 
algorithm to analyze them is derived. This new modeling paradigm should enlarge the class of systems 
amenable at modeling. 

We investigated the influence of both amplitude and autocorrelation time of a extrinsic Sine- Wiener 
noise on: (£) the Michaelis-Menten approximation of noisy enzymatic reactions, which we show to be 
applicable also in co-presence of both intrinsic and extrinsic noise, (ii) a model of enzymatic futile cycle 
and (Hi) a genetic toggle switch. In (ii) and (Hi) we show that the presence of a bounded extrinsic noise 
\^ ' induces qualitative modifications in the probability densities of the involved chemicals, where new modes 

f^i , emerge, thus suggesting the possibile functional role of bounded noises. 

Author Summary 

- ' — \ 

Realistic modeling the dynamics of genetic networks is a non-trivial task that requires choosing a suitable 
level of abstraction. For example, within cells the molecules interacting in a network can be present in 
either small or large numbers. In the former case their time course is typically characterized by wild 
random oscillations closely mimicking the randomness of chemical reactions. In the latter, instead, these 
fluctuations are invisible, due to the "law of large numbers" . So in this case the molecular concentrations 
should theoretically stabilize to nice smooth steady states. 

However, the presence of perturbing external factors may induce noisy fluctuations more or less 
disrupting the theoretical "nice behavior". Surprisingly, this disruption may be constructive: a single- 
function network may gain additional biological functionalities in presence of perturbations. In real world, 
these two kinds of randomness are not separated: proteins of a specific genetic network can be present 
in small number and be perturbed by external noises. This is actually a topic only sparsely explored in 
systems biology. 

A factor that makes complex the study of these phenomena is that external disturbances are classically 
represented through unbounded Gaussian noises, which are actually unrealistic and may induce non- 
biological artifacts. Thus we focus on studying systems with both small number of molecules and external 
bounded noises. 
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After defining an algorithmic framework including the co-presence of intrinsic and extrinsic stochastic 
effects, we apply the developed algorithm to investigate simple motifs of relevance in Systems Biology. 
First, we investigate the possibility of applying the "fast reactions elimination" approximation to en- 
zymatic networks affected by bounded noises. Then, we study a simple model of futile cycle and a 
toggle switch. In both these two cases we show that the addition of extrinsic bounded noises induces the 
emergence of new possible stochastic states. 

Introduction 

Multistability is a key concept of Systems Biology since the first pioneering investigations on the dynamics 
of genetic networks [IH1]. Indeed, it was quite soon understood - both experimentally and theoretically 
- that multiple locally stable equilibria allows for the presence of multiple functionalities, even in small 
groups of interplaying proteins [5Ul4j. 

A second key concept is that the dynamic behavior of a network is never totally deterministic, but it 
exhibits more or less strong stochastic fluctuations due to its interplay with many, and mainly unknown, 
other networks, as well as with various random signals coming from the extracellular world. For long 
time the stochastic effects due these two classes of interactions were interpreted as a disturbance inducing 
undesired jumps between states or, with marginally functional role, as an external initial input directing 
towards one of the possible final states of the network in study. In any case, in the important scenario 
of detcrministically monostablc networks the stochastic behavior under the action of extrinsic noises was 
seen as monomodal. In other words, external stochastic effects were seen similarly as in radiophysics, 
namely as a disturbance more or less obfuscating the real signal, to be controlled by those pathways 
working as a low-pass analog filter [15l[T6]. For these reasons, a number of theoretical and experimental 
investigations focused on the existence of noise-reducing sub- networks [T5J[TT1[IB] - However, it has been 
recently shown the existence of fundamental limits on filtering noise [19) . 

Moreover, if noises were only pure nuisances, there would be an interesting consequence. Indeed, 
in such a case a monostable network in presence of noise should exhibit more or less large fluctuations 
around the unique deterministic equilibrium. In probabilistic languages this means that the probability 
distribution of the total signal (noise plus deterministic signal) should be a sort of "bell" centered more 
or less at the deterministic equilibrium, i.e. the probability distribution should be "unimodal" . However, 
at the end of seventies it became clear in statistical physics that the real stochastic scenario is far 
more complex, and the above-outlined correspondence between deterministic monostability and stochastic 
monomodality in presence of external noise was seriously challenged [20] . Indeed, it was shown that 
many systems that are monostablc in absence of external stochastic noises have, in presence of random 
Gaussian disturbances, multimodal equilibrium probability densities. This counter-intuitive phenomenon 
was termed noise- induced transition [3U] , and it has been shown relevant also in genetic networks [2"Tll2"2"] . 

We above focused mainly on external random perturbations acting on genetic networks. In the 
meantime, experimental studies revealed the other and equally important role of stochastic effects in 
biochemical networks by showing that many important transcription factors, as well as other proteins and 
mRNA, are present in cells with very low concentrations, i.e. with a small number of molecules |23H25j . 
Moreover, it was shown that RNA production is not continuous, but instead it has the characteristics 
of stochastic bursts [26]. Thus, a number of investigations has focused on this internal stochastic effect, 
the "intrinsic noise" as some authors term it |28[l29]. In particular, it was shown - both theoretically 
and experimentally - that also the intrinsic noise may induce multimodality in the discrete probability 
distribution of proteins [22]|30]. However, the fact that intrinsically stochastic systems may exhibit 
behaviors similar to systems affected by extrinsic gaussian noises was very well known in statistical and 
chemical phsyics, where this was theoretically demonstrated by approximating the exact Chemical Master 
Equations with an appropriate Fokker-Planck equation |31H33| . an approach leading to the Chemical 
Langcvin Equation [34] . 
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Thus, after that for some time noise was mostly seen as a nuisance, more recently it has finally been 
appreciated that the above-mentioned and other noise-related phenomena may in many cases have a 
constructive, functional role (e.g. see |35U36j and references therein). Indeed, for example, noise-induced 
multimodality allows a transcription network for reaching states that would not be accessible if the noise 
was absent [22[|35[|36] . Phenotype variability in cellular populations is probably the most important 
macroscopic effect of intracellular noise- induced multimodality |35j . 

In Systems Biology, from the modeling point of view Swain and coworkers |24| were among the 
first to study the co-presence of both intrinsic and extrinsic randomness, by stressing the synergic role in 
modifying the velocity and average in the context of the basic network for the production and consumption 
of a single protein, in absence of feedbacks. These and other important effects were shown, although 
nonlinear phenomena such as multimodality were absent. The above study is also remarkable since: 
(i) it has stressed the role of the autocorrelation time of the external noise and, differently from other 
investigations, (ii) it has stressed that modeling the external noise by means of a Gaussian noise, either 
white or colored, may induce artifacts. In fact, since the perturbed parameters may become negative, the 
authors employed a lognormal positive noise to model the extrinsic perturbations. In particular, in |24j 
a noise obtained by exponentiating the classical Orenstin-Uhlenbeck noise was used [20) . 

From the data analysis point of view, You and collaborators [37] and Hilfingcr and Paulsson [38] 
recently proposed interesting methodologies to infer by convolution the contributions of extrinsic noise 
also in some nonlinear networks, including a synthetic toggle switch [37] . 

Our general aim here is manifold. Indeed, we want to start by investigating the co-presence of both 
extrinsic and intrinsic randomness in nonlinear genetic networks in the important case where the external 
perturbations are not only non-Gaussian, but also bounded. Indeed, by imposing a bounded extrinsic 
noise we increase the degree of realism of a model, since the external perturbations must not only preserve 
the positiveness of reaction rates, but must also be bounded. Moreover, it has also been shown in other 
contexts such as oncology and statistical physics that: (i) bounded noises deeply impact on the transitions 
from unimodal to multimodal probability distribution of state variables [39H43] and (ii) the dynamics of 
a system under bounded noise may be substantially different from the one of systems perturbed by other 
kinds of noises, for example there is dependence of the behavior on the initial conditions [40 ] I42 | . 

This approach opens a number of questions, two of which we shall try to assess here. The first question 
is to identify a suitable mathematical framework to represent mass-action biochemical networks perturbed 
by bounded noises (or simply left-bounded), which in turn can depend on the state of the system. To this 
extent, in the first part of this work we derive a master equation for these kinds of systems in terms of the 
differential Chapman-Kolgomorov equation (DCKE) [31J44J and propose a combination of the Gillespie's 
Stochastic Simulation Algorithm (SSA) [27l[28] with a state-dependent Langevin system, affecting the 
model jump rates, to simulate these systems. 

The second question relates to the possibility of extending, in this "doubly stochastic" context, the 
Michaclis-Mcnten Quasi Steady State approximation (QSSA) for enzymatic reactions [45] ■ We face the 
validity of the QSSA in coprcsence of both types of noise in the second part of this work, where we 
numerically investigate the classical Enzymc-Substratc-Product network. The application of QSSA in 
this network has been recently investigated by Gillespie and coworkers when extrinsic noise is absent [46] . 
Based on our results, we propose the extension of the above structure also to more general networks than 
those ruled by the rigourous mass-action law via a stochastic QSSA. 

Finally, in the third part we investigate the interplay between intrinsic randomness and extrinsic 
bounded noises in two cases of interest in biology: (i) a futile cycle [22] and (ii) a genetic toggle switch [7], 
which is a fundamental motif for cellular differentiation and for other switching functions. As expected, 
the co-presence of both intrinsic stochasticity and bounded extrinsic random perturbations suggests the 
presence of possibly unknown functional roles for noise in both networks. The described noise-induced 
phenomena are shown to be strongly related to physical characteristics of the extrinsic noise such as the 
noise amplitude and its autocorrelation time. 
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Table 1. Analytical form of the propensity functions [27] . 



Methods 

Noise-free stochastic chemically reacting systems 

We start by recalling the Chemical Master Equation and the Stochastic Simulation Algorithm (SSA) by 
Doob and Gillespie J27J[25]. Systems where the jump rates are time-constant are hereby referred to as 
stochastic noise-free systems. We consider a well stirred system of molecules belonging to N chemical 
species {Si, ... , Sjv} interacting through M chemical reactions Ri, . . . , Rm- We represent the (discrete) 
state of the target system with a A-dimcnsional integer- valued vector X(i) = (Xi(t), . . . , .Xjv(t)) where 
Xi(t) is the number of molecules of species Si at time t. To each reaction Rj is associated its stoichiometric 
vector Vj = ■ • ■ , v N,j), where Vij is the change in the Xi due to one Rj reaction. The stoichiometric 
vectors form the N x M stoichiomctry matrix D = \v\ ... Vm\ ■ Thus, given X(t) = x the firing 

of reaction Rj yields the new state x + Vj . A propensity function dj (x) [27J[2H] is associated to each Rj 
so that a,j(jx.)dt, given X(i) = x, is the probability of reaction Rj to fire in state x in the infinitesimal 
interval [t, t + dt). Table [1] summarizes the analytical form of such functions [27] . For more generic form 
of the propensity functions (e.g. Michaelis-Menten, Hill kinetics) we refer to |49j . 

We recall the definition of the Chemical Master Equation (CME) [27 l l28 ] l47 ] l48 ] describing the time- 
evolution of the probability of a system to occupy each one of a set of states. We study the time-evolution 
of X(i). assuming that the system was initially in some state xo at time to, i.e. X(to) = xq. Wc denote 
with V(x,t | xo,to) = "P(x,i | uj) the probability that, given X(to) = xo, at time t it is X(i) = x. From 
the usual hypothesis that at most one reaction fires in the infinitesimal interval [t, t + dt), it follows that 
the time-evolution of 'P(x, t \ lo) is given by the following partial differential equation termed "master 
equation" 

M 

dtP(x,t \u)=^2v{x-Vj,t\ w) aj -(x -Vj) - V(x,t \ w)a,-(x) . (1) 
i=i 

The CME is a special case of the more general Kolmogorov Equations [50], i.e. the differential equations 
corresponding to the time-evolution of stochastic Markov jump processes. As it is well known, the 
CME can be solved analytcally only for a very few simple systems, and normalization techniques are 
sometimes adopted to provide approximate solutions [ST] . However, exact algorithmic realization of the 
process associated to the CME are possible by using the Doob-Gillcspie Stochastic Simulation Algorithm 
(SSA) [27l[28l[47l[48] , summarized as Algorithm [TJ Although equivalent formulations exist [27J[28l[52] , as 
well as some approximations [H][S31[S5 , here we consider its Direct Method formulation without loss of 
generality. 

The SSA is an exact dynamic Monte-Carlo method describing a statistically correct trajectory of a 
discrete non-linear Markov process, whose probability density function is the solution of equation ([TJ [55] . 
The SSA computes a single realization of the process X(t), starting from state xo at time to and up to 
time T. Given X(t) = x the putative time r for the next reaction to fire is chosen by sampling an 
exponentially distributed random variable, i.e. r ~ £xp(ao(x)) where ao(x) = J2j=i a j'( x ) an d ~ denotes 
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Algorithm 1 SSA (t , x , T) 

1: set x xo and t to; 
2: while t <T do 

3: o (x) <- Ejii°j( x ); 

4: let n, r 2 ~ C/[0,1]; 
5: r <- ^(x)- 1 !!!^; 1 ); 

6: let j such that J2l=i a *( x ) < r 2 • ao( x ) < X)<=i a i( x ); 
7: set x <— x + Vj and t <— t + r; 

8: end while 



the equality in law between random variables. The reaction to fire Rj is chosen with weighted probability 
a,j (x) / ao (x) , and the system state is updated accordingly. 

The correctness of the SSA comes from the relation between the jump process and the CME [27l[55] . 
In fact, the probability, given X(£) = x, that the next reaction in the system occurs in the infinitesimal 
time interval [t + r, t + t + dr), denoted V(t | x, t), follows 



V(t |M)=X) P(t, j I x, t) = ao (x) exp {- £ a (t')dt'^j = a (x) e - ao(x)T 



(2) 



since V(t, j x, t) = Oj(x) exp (— ao(x)r) is the probability distribution of the putative time for the next 
firing of Rj, and the formula follows by the independency of the reaction firings. Notice that in equation 
@ &o(t') represents the propensity functions evaluated in the system state at time t' > t, i.e. as if they 
were time-dependent functions. In the case of noise-free systems that term evaluates as cto(x) for any 
t G \t,t + t], i.e. it is indeed time-homogenous whereas in more general cases it may not, as we shall 
discuss later. Finally, the probability of the reaction to fire at t + r to be Rj follows by conditioning on 
j, that is 

v( - , .s = 'PjrJ | x,f) = « 3 (x) exp(-a (x)r) = aj -(x) 
U|T '' J V{r\*,t) ao(x)exp(-a (x)T) a (x) ' W 

Noisy stochastic chemically reacting systems 

We now introduce a theory of stochastic chemically reacting systems with (un)bounded noise in the jump 
rates by combining Stochastic Differential Equations and the SSA. Here we consider a system where each 
propensity function may be affected by a extrinsic noise term. In general, such a term can be either a 
time or state-dependent function, and the propensity function for reaction Rj reads now as 

a j (x,t)=a j (x)L*(t), (4) 

where Oj (x) is a propensity function of a type listed in Table [TJ The noisy perturbation term L* (t) is 
positive and bounded by some Cj < +oo, i.e. 

< L*(t) < Cj (5) 

so we are actually considering both bounded and right-unbounded noises, i.e. Cj = +oo. In the former 
case we say that the j-th extrinsic noise is bounded, in the latter that it is left-bounded. Note that in 
applications we shall mainly consider unitary mean perturbations, that is 

(L*(t)) = l. 

We consider here that the extrinsic noisy disturbance L* (t) is a function of a more generic E-dimensional 
noise with 1 < S < M so we write L*{t) = Lj(£(t)) and equation (@| reads as 

a,(x,i) =a,(x)L,(£(t)). (6) 
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Notice that the use of a vector in equation provides the important case of multiple reactions sharing 
the same noise term. i.e. the reactions may be affected in the same way by a unique noise source. In 
equation © Lj is a continuous functions Lj : IR S — > K and G M s is a colored and, in general, 
non-gaussian noise that may depend on the state X(f) of the chemical system. The dynamics of £(t) is 
described by a E-dimensional Langevin system 

Z'(t)=m,X(t))+g(Z,X(t))ri(t). (7) 

Here, r\ is a S-dimensional vector of uncorrclated white noises of unitary intensities, g is a S x X matrix 
which we shall mainly consider the be diagonal and /, gu,k '■ ^ S x R w — ► R S - When does not directly 
depend on X(i), i.e. the extrinsic noise depends on an external source, which is the kind of noise we 
mainly consider, equation Q reduces to 

e(t) = f(0 + g(Ov(t)- (8) 

We stress that the "complete" Langevin system in equation ([7]) is not a mere analytical exercise, but 
it has the aim of phenomenologically modeling extrinsic noises that are not totally independent of the 
process in study. 



The Chapman-Kolmogorov Forward Equation 

When a discrete-state jump process as one of those described in previous section is linked with a continuous 
noise the state of the stochastic process is the vector 

z = (x, f) where x G N N , £ G K s , (9) 

and the state space of the process is now x K s . Our total process can be considered MS (1 particular 
case of the general Markov process where diffusion, drift and discrete finite jumps are all co-present for all 
state variables [3T|l44] . For this very general family of stochastic processes the dynamics of the probability 
of being in some state z at time t, given an initial state Zq at time to shortly denoted as w, is described 
by the differential Chapman-Kolgomorov equation (DCKE) |31U44j . whose generic form is 

d t v(z,t | w) =-J2^ J A j (z,t)v(z,t | w) +^Y / d ZuZj B id (z,t)v(z,t | w) (10) 

j i,j 

+ J \\¥{z | h, t)v(z, t | w) - W{h | z, t)v(h, t \ uj^ 

Here Aj forms the drift vector for z, Bij(z,t) the diffusion matrix and W the jump probability. For an 
elegant derivation of the DCKE from the integral Chapman-Kolgomorov equation [50] we refer to |44| . 
This equation describes various systems, in fact we remind that (i) the Fokker-Planck equation is a 
particular case of the DCKE without jumps (i.e. W(z \ h, t) = 0), (ii) the CME in equation ((1) is the 
DCKE without brownian motion and drift (i.e. A(z,t) = and B(z,t) = 0), (Hi) the Liouville equation 
is the DCKE without brownian motion and jumps (i.e. A(z,t) = and W(z \ h, t) = 0) and (iv) the 
ODE with jumps correspond to the case where only diffusion is absent (i.e. B(z,t) = 0). 

We stress that, at the best of our knowledge, this is the first time where a master equation for 
stochastic chemically reacting systems combined with bounded noises is considered. Let 

p((x,C),* I (xo, &),*(>) =V(m,t\u) (11) 

be the probability that at time t it is X(t) = x and = £, given X(£o) = xo and £(io) = Co- The 
time-evolution of V(z, t \ u>) is equation (|10j) where drift and diffusion are given by the Langevin equation 
0, that is 

A = /(£,x) S = .9 T x.g (12) 
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with x the standard vector multiplication and g T the transpose of g. Moreover, since only finite jumps 
are possible, then the jump functions and diffusion satisfy 

d XiXj B itj (z,t) = W((x,0 I (x,r),i) =0 (13) 

for any i,j = 1,...,N, and noise £* £ M s . Summarizing, for the systems we consider the DCKE in 
equation (fTQ|) reads as 

M M M 

d t p((x,0,t I w) =-^9 z ,/ J (^x)7'((x ) e),i I w) +2EE^ fl yK' 3t ) ? '(( x '0 1 i I w) (14) 

3=1 <=1 j=l 

A/ A/ 

This equation is the natural generalization of the CME in equation (fTJ), and completely characterize noisy 
systems. As such, however, its realization can be prohibitively difficult and is hence convenient to define 
algorithms to perform the simulation of noisy systems. 



The SSA with Bounded Noise 

We now define the Stochastic Simulation Algorithm with Langevin Noise (SSAn). The algorithm performs 
a realization of the stochastic process underlying the system where a (generic) realization of the noise is 
assumed. As for the CME and the SSA, this corresponds to computing a realization of a process satisfying 
equation (Ti"4|) . The SSAn takes inspiration from the (generic) SSA with time-dependent propensity 
functions [56] as well as the SSA for hybrid deterministic/stochastic systems [57Tl60j . thus generalizing 
the jump equation ((2|) to a time inhomogeneous distribution, which we discuss in the following. 
For a system with M reactions the time evolution equation for X(i) is 

M 

dX(t) (15) 

where {Nj(t) \ t > to} is the stochastic process counting the number of times that Rj occurs in [to, t] with 
initial condition Nj(to) = 0. For Markov processes Nj(t) is an inhomogeneous Poisson process satisfying 

v(Nj(t + dt) - Nj(t) = 1 [ x) = a,j(x.,t)dt = aj(x)Lj(£(t))dt (16) 

when X(t) = x. In hybrid systems this is is a doubly stochastic Poisson process with timc-dependent 
intensity, in our case this is a Cox process [61][62] since the intensity itself is a stochastic process, 
i.e. it depends on the stochastic noise. More simply, in noise-free systems, this equation evaluates as 
aj(x)dt, thus denoting a time homogeneous Poisson process. As in [57,59,60,63,641 such a process ca be 
transformed in a time homogenous Poisson process with parameter 1, and a simulation algorithm can be 
exploited. Let us denote with Tj(t) the time at next occurrence of reaction Rj after time t, then 

v(Tj(t) e[t,t + dt] | x) = aj(x)i 3 -(^(*))di + o(dt) (17) 

follows by equation (|16[) and higher order terms vanish by the usual hypothesis that the reaction firings 
are locally independent, as in the derivation of equation (fTJ. Given the system to be in state x at time 
t, the transformation 



rt+r 

A j (t,t + T)=a j (x) J L,(at'))dt' 



(18) 
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Algorithm 2 SSAn (t , x , T) 

1: set x xo and t <— to! 
2: while t < T do 

3: let n, r 2 ~ U[0,1]; 

4: find r by solving J2f=i A j(t + T ,t) = M^ 1 ) an d for £ [t,t + r]; 

5: let j such that a *( x ' * + T ) < r 2 Sfeli a fc( x - * + r ) < Z)i=i a *( x ' f + T )? 

6: set x x + Vj and i <— i + r; 
7: end while 



which is a monotonic (increasing) function of t is used to determine the putative time for Rj to fire. 
Given a sequence Tj k of independent exponential random variables with mean 1 for j = 1, . . . , M and 
fceN, equation (fT6)l implies that 

oo 

^•t^E^^w}' ( 19 ) 

This provides that, if the systems is in state X(i) = x, then the next time for the next reaction firing of 
Rj is the smallest time r > such that 

Aj(t + T,t) = r (20) 

with r ~ fxp(l), and thus the next jump of the overall system is taken as the minimum among all 
possible times, that is by solving equality 

M m t+T 

E^(* + M)=5>W / L Mt'))dt'=r (21) 

with r ~ £xp(l). This holds because minjTj | j = 1,...,M} is still exponential with parameter 
^2jLi Aj(t + t, t) and the jumps are independent. We remark that for a noise-free reaction Aj(t + t, t) = 
ra 3 -(x), thus suggesting that the combination of noisy and noise- free reactions is straightforward. The 
index of the reaction to fire is instead a random variable following 

, , s a i (x)L i (e(t + T)) , , 

V(j r; x, t) = ./ V ; JV ^ V ^— . 22 

+ r)) 

The SSAn is Algorithm^ Step 4 is the (parallel) solution of both equation ([2Tj) and Langcvin system 

0, step 5 samples values for j according to equation (f22|) . As far as step 4 is concerned, it is worth 
nothing that given X(i) = x for any r the Langevin equation ([7} depends only on £(t) and the constant 
x. To this extent, a single trajectory of the vectorial noise process in [t,t + t] is 

£t,r = {(*, £(*))} U {(* s , £(* s )) I * < t" < t + t} U {(t + r, £(i + r))} . (23) 

This is a discretization of a continuous noise, thus inducing an approximation, but is in general the only 
possible approach. To reduce approximation errors the maximum size of the jump in the noise realization, 

1. e. the noise granularity A s = i s+1 — t s , should be much smaller than the minimum autocorrelation time 
of the perturbing stochastic processes Lj(£(t)). 

Finally, the integral in equation (|21j) . evaluated in step 4, is a conventional Lebesguc integral since 
the perturbation Lj(£(t)) is a colored stochastic process ^53- As an example, given £ t T by a linear 
interpolation scheme it holds 

LMs))ds* J2 A B (^{L%Lp- 1 } + ±\L' j -Lf 1 \\ (24) 
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where = Lj(^(t s )) for t s S £ t T and A s the noise granularity. 
Extension to non mass-action nonlinear kinetic laws 

Large networks with large chemical concentrations, i.e. characterized by deterministic behaviors, are 
amenable to significant simplifications by means of the well known Quasi Steady State Approximation 
(QSSA) [SISHlSSlinB] ■ The validity conditions underlying these assumptions are very well-known in the 
context of deterministic models |45j , despite not much being known for the corresponding stochastic 
models. Recently, Gillespie and coworkers [46] showed that, in the classical Michaelis-Menten Enzyme- 
Substratc-Product network, a kind of Stochastic QSSA (SQSSA) may be applied as well, and that in 
such its limitations are identical to the deterministic QSSA. Thus, it is of interest to consider SQSSAs 
also in our "doubly stochastic" setting, even though possible pitfalls may arise due to the presence of 
the extrinsic noises. As an example, in Section „ we will present the results of the numerical experiments 
similar to those of [46], with the purpose of validating the SQSSA for noisy Michaelis-Menten enzymatic 
reactions. 

Of course, in a SQSSA not only the propensities may be nonlinear function of state variables, but 
they may depend nonlinearly also on the perturbations, so that instead of the elementary perturbed 
propensities we shall have generalized perturbed propensities of the form 

where ip is a vector with elements ijij = Lj(£) for j = 1, . . . , M. This makes possible, within the above 
outlined limitation for the applicability of the SQSSA, to write a DCKE for these systems as 

M 

E 

i=i 

M M 



$7>((x,0,tk) w x)7>((x,0,t I <r) (25) 
3=1 
M M 

+ 2EE 9 ^ s « & ^ V (( x - ' ^ I a ) 

i=l 0=1 
M 

+ E V (( x ~ v i> 0. * I <*J ( x - V 3 , V- (*)) 



3=1 

M 
3=1 

As far as the simulation algorithm is concerned, it remains quite close to Algorithm [5] provided that the 
jump times are sampled according to the following distribution 

p(r I x, t) = ctj(x, tp(t + r)) exp I — / aj(x,t(}(k))dk\ . (26) 



Results 

We performed SSAN-based analysis of some simple biological networks, actually present in most complex 
realistic networks. We start by studying the legitimacy of the stochastic Michaelis-Menten approximation 
of when noise affects enzyme kinetics |46j . Then we study the role of the copresence of intrinsic and etrinsic 
bounded noises in a in a model of enzymatic futile cycle [22] and, finally, in a bistable "toggle switch" 
model of gene expression [13j[73] . All the simulations have been performed by a Java implementation of 
the SSAn (freely downloadable at http://sites.google.com/site/giuliocaravagna/) running on a 



cluster of 15 dual-core nodes with 2.0 Ghz processor and 1 GB of memory. 



10 



The Sine- Wiener noise |39j. The bounded noise /j,(t) that we use in our simulations is obtained by 
applying a bounded continuous function h : ffi — > R to a random walk W, i.e. W'(t) = rj(t) with rj(t) a 
white noiscQ. We have 

m = h(w) 

so that for some /? £ R it holds — (3 < h(W) < (3. The effect of the truncation of the tails induced by the 
approach here illustrated is that, due to this "compression", the stationary probability densities of this 
class of processes satisfy 

V(ji = \p\) = +co. 

Probably the best studied bounded stochastic process obtained by using this approach is the so-called 
Sine- Wiener noise [39], that is 

fi(t) = /3 sin (^W(t)j (27) 

where (3 is the noise intensity and r is the autocorrelation time. The average and the variance of this 
noise are: 

<M*)> = W) 2 ) = l3 2 /2 

and its autocorrelation is such that IMI 



(3 2 ( -z 
{n{t)n(t + z )) = — cxp ( — 



1 — cxp 



-At 

T 



Note that, since we mean to use noises of the form 1 4- /i(t), i.e. the unitary- mean perturbations in 
equation ©, then the noise amplitude must be such that < (3 < 1. 
For this noise, the probability density is the following [77] : 

TO 



By these properties, this noise can be considered a realistic extension of the well-known symmetric 
dichotomous Markov noise [17] a(t), whose stationary density is \ (5(a — f3) + S(a + (3)), for a e R and 5 
the Dirac delta function. 



Enzyme kinetics 

Enzyme-catalyzed reactions are fundamental for life, and in deterministic chemical kinetics theories are 
often conveniently represented in an approximated non mass-action form, the well-known Michaclis- 
Mcntcn kinetics [311351135] ■ Such approximation of the exact mass-action model is based on a Quasi 
Steady-State Assumption (QSSA) [3SE1], valid under some well known conditions. In [32] it is studied the 
legitimacy of the Michaclis-Mcnten approximation of the Enzyme-Substrate-Product stochastic reaction 
kinetics. Most important, it is shown that such a stochastic approximation, i.e. the SQSSA previously 
discussed, obeys the same validity conditions for the deterministic regime. This suggests the legitimacy 
of using - in case of low number of molecules - the Gillespie algorithm not only for simulationg mass- 
action law kinetics, but more in general to simulate more complex rate laws, once a simple conversion of 

lr The underlying white-noise process W(t) is generated at times {ti | i > 0} according to the recursive schema W(t;+i) = 
W(ti) + Ti \/Aw with initial condition W(to) = tq. Here n £ A^(0, 1) and Aw = ti+l — ti for i > is its discretization 
step; it has to satisfy Aw <C r and thus we chose to be Ayy i=s r/100. Notice that, as intuitive the noise autocorrelation is 
expected to deeply impact on the simulation times. 
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-1 1 1 

-1 1 

1 -1 -1 

1 



ai = ciE ■ S 

0,2 = C2ES 

03 = csES 



-1 


V 1 



ai = 



Vmax $ 

K M + S 



Table 2. Exact model (left) and Michaelis-Menten approximation (right) of enzymatic reactions: the 
stoichiometry matrixes (rows in order E, S, ES, P) and the propensity functions. 



(i) (E, S, ES, P) = (1, 10, 0, 0) (ii) (E, S, ES, P) = (100, 10, 0, 0) 




10 20 30 40 50 1 2 3 4 5 



Figure 1. Noise-free Enzyme-Substrate-Product system. Averages of 1000 simulations for P, 
plotted with its standard deviation (dotted), for both exact and approximated Michaelis-Menten 
models. We have set c± = C3 = 1 and c 2 = 10 and the initial configuration is (i) in left 
(E, S, ES, P) = (1, 10, 0, 0) and (ii) in right (E, S, ES, P) = (100, 10, 0, 0). 

deterministic Michaelis-Menten models is performed and provided - of course - that the SQSSA validity 
conditions arc fulfilled. 

In this section we investigate numerically whether the Michaelis-Menten approximations and the 
stochastic results obtained in [35] still hold true in case that a bounded stochastic noise perturb the 
kinetic constants of the propensities of the exact mass-action law system Enzyme-Substrate-Product. 
Let E be an enzyme, S a substrate and P a product, the exact mass-action model of enzymatic reactions 
comprises the following three reactions 

E + S^ES ES^E + S ES^E + P 

where c\, Ci and C3 are the kinetic constants. The network describes the transformation of substrate S 
into product P , as driven by the formation of the enzyme-substrate complex ES, which is reversible. 
The deterministic version of such reactions is 

S' = -c x S ■ E + c 2 ES E' = - Cl S ■ E + (c 2 + c 3 )ES (28) 

ES' = aS-E- (c 2 + c 3 )ES P' = c 3 ES, 

where we write 5 • E to distinguish the multiplication of E and S from complex ES. By the relations 

E T = E(t)+ES(t) P(t) = P(t ) + S(t )-(S(t) + ES(t)) (29) 

a QSSA reduces to one the number of involved equations. Indeed, since ES is in quasi-steady-state, i.e. 
ES' = 0, then 

pl ^ V max where V max = c 3 E T and K M = C% + ° 3 . (30) 

K M + S ci 

Here Km is termed the Michaelis-Menten constant. In practice, the QSSA permits to reduce the three- 
reactions model to the single-reaction model 

S^P 
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(up) n = 1, ft = 0.5 (down) t» = 1, ft = 1 (up) n = 1, ft = 0.5 (down) r< = 1, ft = 1 




10 20 30 40 50 1 2 3 4 5 



(up) n = 5, ft = 0.5 (down) r 4 = 5, ft = 1 (up) r t = 5, ft = 0.5 (down) r< = 5, ft = 1 




10 20 30 40 50 1 2 3 4 5 



Figure 2. Stochastically perturbed Enzyme-Substrate-Product system. Averages of 1000 
simulations for P, plotted with its standard deviation (dotted), for both exact and approximated 
Michaelis-Menten models. Parameters are as in Figure [T](z) in the top four panels, and as in Figured] 
(ii) in the bottom four. Independent Sine- Wiener noises are present in all the reactions. For i = 1,2, 3, 
in the left column panels Tj = 1, in the right r, = 5, in top = 0.5 and in bottom Pi = 1. 

with non mass-action non linear rate (V max S)/(KM + S). In [55] the condition 

E T « So + K M (31) 

is used to determine a region of the parameters space guaranteeing the legitimacy of the Michaelis-Menten 
approximation. When condition (|31|) holds, a separation exists between the fast pre-steady-state and the 
slower steady-state timescales [66] and the solution of the Michaelis-Menten approximation closely tracks 
the solution of the exact model on the slow timescale. 

Here we show that the same condition is sufficient to legitimate the Michaelis-Menten approximation 
with bounded noises arbitrarily applied to any of the involved reactions. We start by recalling the result 
in [46] about the noise-free models given in Table [2] We considered two initial conditions: (i) one with 
10 copies of substrate, 1 enzyme and complexes and products, and (ii) one with 10 copies of substrate, 
100 enzyme and complexes and products. As in j46| we set c\ = C3 = 1 and C2 = 10; notice that the 
parameters are dimensionless and, more important, in (i) they satisfy condition (|3ip since Et = 1 and 
So + Km = 21, in (ii) no. In Figure [T] we reproduced the results in |46] for (i) in right panel and (ii) in 
left. As expected, in (i) the approximation is valid on the slow time-scale, and not valid in the fast, i.e. 
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(up) ti = 1, /3i — 0.5 (down) n — 1, /3i — 1 (up) T2 — 1, /?2 — 0.5 (down) T2 = 1, /?2 = 1 





(up) r 3 = 1, ^3 = 0.5 (down) t 3 = 1, 3 = 1 (up) n = 5, /3i = 0.5 (down) n = 5, 01 = 1 




(up) t 2 = 5, £2 = 0.5 (down) r 2 = 5, /3 2 = 1 (up) t 3 = 5, /3 3 = 0.5 (down) r 3 = 5, /3 3 = 1 





Figure 3. Stochastically perturbed Enzyme-Substrate-Product system. Averages of 1000 
simulations for P with dotted standard deviation, for both exact and approximated Michaelis-Menten 
models with parameters as in Figure Q] (i). Here single noises with two intensities and different 
autocorrelations are used. The non-zero parameters are reported in top captions. 



for t < 3, in (ii) it is not valid also in the slow time-scale. 

If noises are considered the models in Table [5] change accordingly. So, for instance when independent 
Sine- Wiener noises are applied to each reaction, the exact model becomes 



ai = ci 



1 + ft sin ^v^7^PFi(i)jJ E-S a 2 = c 2 [l + ft sin (^2/7 2 W 2 (t) 

a 3 = c 3 [l + ft sin (v/27^ 3 (t))] ES 
and the Michaelis-Menten constant becomes the time-dependent function 



ES 



K* M (t) 



C2 


1 + ft sin 


[^2/r 2 W 2 {t)) 


+ c 3 


1 + ft sin 




f'l 


1 + ft sin 


^2/r 1 W 1 (t) J 







Notice that the nonlinear approximated propensity ai(t) is now time-dependent, and, moreover, it de- 
pends nonlincarly on the noises affecting the system. 
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(up) Ti = 1, 01 = 0.5 (down) ti = 1, 0i = 1 (up) Ti = 5, /3i = 0.5 (down) ri = 5, 01 = 1 




Figure 4. Stochastically perturbed Enzyme-Substrate-Product system. Averages of 1000 
simulations for P, plotted with its standard deviation (dotted) in the fast time-scale < t < 5 for both 
exact and approximated Michaelis-Menten models with parameters as in Figure Q] (i). Here we use 
independent Sine- Wiener noises with n = 1 in left and t\ = 5 in right panels. For every parameters 
configuration both low, i.e. 0.5, and high, i.e. 1, noise level intensities are used. 

Thus condition (|31|) becomes time-dependent and we rephrase it to be 

E T <C S + (K* M (t)) . (32) 

Note that if ft > then {K^{t)) ^ K M , whereas if ft = then (K* M (t)) = K M - 

Each of the shown figures is the result of 1000 simulations for model configuration where the simulation 
times, which span from few seconds to few minutes, depend on the noise correlation. When the same 
system of Figure [T] (i) is extended with these noises the approximation is still valid, as shown in the 
top panels of Figure [2j In addition, the approximation is not valid when condition (|32|) does not hold, 
as shown in the bottom panels of Figure [21 as it was in Figure Q] (ii). Notice that in there we use two 
different noise correlations, i.e. Ti = 1 in the left and Tj = 5 for i = 1,2,3 in the right column panels, 
thus mimicking noise sources with quite different charateristic kinetics. Also, we set two different noise 
intensities, i.e. ft = 0.5 in top panels and ft = 1 (maximum intensity) in bottom panels, whereas all 
the other parameters are as in Figure [U Summarizing, we get a complete agreement between enzymatic 
reactions with/without noise, independently on the noise characteristics when it affects all of the reactions. 

To strengthen this conclusion it becomes important to investigate whether it still holds when noises 
affects only a portion of the network and, also, whether it holds on the fast time-scale. 

As far as the number of noises is concerned, we investigated various single-noise configurations in 
Figure [3] In there we used a single noise, i.e. two out of the three noises have intensity, with both low 
and high intensities, i.e. 0.5 and 1. Also, in that figure we vary the noise correlation time as r G [1,5]. 
As hoped, the simulations show that the approximation is legitimate in the slow time-scale for all the 
various parameter configurations, thus independently on the presence of single or multiple noises. 

Finally, as far as the legitimacy of the approximation in the fast time-scale is concerned, i.e. t G [0, 5], 
our simulations show a result of interest: if the noise correlation is small compared to the reference fast 
time-scale and if single noises are considered the noisy Michaelis-Menten approximation performs well 
also on the fast time-scale. We remark that this was not the case for the analogous noise-free scenario in 
Figure [T] (i). In support of this we plot in Figure [3] the fast time-scale for rj = 1 and r, = 5 for the single 
noise model with a noise in the enzyme-substrate complex formation, i.e. ft = ft = 0. Similar evidences 
were found in the configurations plotted in Figure [3] (not shown) . 
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kiE - So a,2 = k-iESa 
= k 2 ES a4 = k 3 FSi 
k- 3 FSi a 6 = k±FSi 



Table 3. Futile cycle model. The noise- free enzymatic futile cycle [22]: the stoichiometry matrix 
(rows in order So, E, ESq, Si, F, FSi) and the propensity functions. 

Futile cycles 

In this section we consider a model of futile cycle, as the one computationally studied in [22) . The model 
consists of the following mass-action reactions 

E + S Q ^ ES Q ESo ^> E + S Q ESo^E + Si 

F + Si ^ FSi FS*i 3-^ F + Si FSi^F + So 

where E and F are enzymes, So and Si substrate molecules, and ESq and FSi the complexes enzyme- 
substrate. Futile cycles are an unbiquitous class of biochemical reactions, acing as a motif in many signal 
transduction pathways [68] . 

Experimental evidences related the presence of enzymatic cycles with bimodalities in stochastic chem- 
ical activities [55]. As already seen in the previous section, Michaelis-Menten kinetics is not sufficient to 
describe such complex behaviors, and further enzymatic processes are often introduced to induce more 
complex behaviors. For instance, in deterministic models of enzymatic reactions feedbacks are necessary 
to induce bifurcations and oscillations. Instead, in [22] it is shown that, although the determinstic version 
of the model has a unique and attractive equilibrium state, stochastic fluctuations in the total number 
of E molecules may induce a transition from a unimodal to a bimodal behaviour of the chemicals. This 
phenomenon was shown both by the analytical study of a continuous SDE model where the random 
fluctuations in the total number of enzyme E (both free and as a complex with S) is modeled by means 
of a white gaussian noise on the one hand, and in a totally stochastic setting on the other hand. In the 
latter case it was assumed the presence of a third molecule N interacting with enzyme E according to 
the following reactions 

N + N ^ E + N E + N N + N 

N^E E^N. 

By using N the stochastic model results to be both quantitatively and qualitatively different from the 
deterministic equivalent. These differences serve to confer additional functional modalities on the enzy- 
matic futile cycle mechanism that include stochastic amplification and signaling, the characteristics of 
which depend on the noise. 

Our aim here is to investigate whether bounded noises affecting the kinetic constant, and thus not 
modifying the topology of the futile cycle network, may as well induce transition to bimodality in the 
system behavior. To this aim, here we analyze three model configurations: (i) the noise-free futile cycle, 
namely only the first six reactions, (ii) the futile cycle with the external noise as given by N and (iii) 
the futile cycle with a bounded noise on the binding of E and So, i-e. the formation of ESo, an d N is 
absent. 



16 



noise- free 



futile cycle with ./V 




cycle with bounded noise with r — 0.1 
(top) = 0.5 (bottom) = 1 



cycle with bounded noise with r — 1 
(top) P = 0.5 (bottom) P = 1 
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Figure 5. Stochastic models of futile cycles. Single run and averages of 1000 simulations for 
substrate So of the futile cycle models. In top panels the noise-free (left) and the cycle unimodal noise 
as N (right). In bottom panels the cycle with bounded noise and r = 0.1 in (left) and t = 1 in (right). 
In both cases in the top panel the noise intensity is f3 = 0.5 (top) and j3 = 1 (bottom). The initial 
configuration is always (So, E, ES , Si,F, FSi,N) = (0, 20, 0, 2000, 50, 0, 10) and the kinetic parameters 
are k\ = 40, fc_i = k 2 = 10000, k 3 = 200, k_ 3 = 100, k 4 = 5000 for the noise-free and the bounded noise 
case, and k 5 = fc 6 = 10, /c_ 5 = 5 and fc_ 6 = 0.2 [2"2"] . 



In Table [3] the noise-free futile cycle is given as a stoichiomctry matrix and 6 mass-action reactions. 
The model simulated in |22| is obtained by extending the model in the table with a stoichiometry matrix 
containing N and four more mass-action reactions. For the sake of shortening the presentation we omit 
to show them here. The model with a bounded noise in a\ is obtained by defining 



ai(t) = k\E ■ Sq 



1 + (3 sin 




We simulated the above three models according to the initial condition used in [52] (5*0, E, ES , Si,F, FSi) 
(0, 20, 0, 2000, 50, 0) which is extended to account for 10 initial molecules of N, when necessary. The ki- 
netic parameters are dimensionless and defined as k\ = 40, = k 2 — 10000, k 3 = 200, fc„3 = 100, 
fc 4 = 5000 for the noise-free and the bounded noise case, and k$ = fc 6 = 10, fc_ 5 = 5 and fc_6 = 0.2 when 
the unimodal noise is considered [22] . Furthermore, when the bounded noise is considered the autocor- 
relation is chosen as r G [k^ 1 , 1] = [0.1, 1] according to the highest rate of the reactions generating the 
unimodal noise. 

In Figure [5] a single run and averages of 1000 simulations for the futile cycle models are shown. In 
this case the simulation times span from 20 -f- 30 rain to 60 80 min, thus making the choice of good 
parameters more crucial than in the other cases. In Figure [S] the substrate 5*0 is plotted, and 5*1 behaves 
complementarily. In top panels the noise-free (top) and the cycle unimodal noise as N (bottom). In 
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noisc-frce cycle cycle with species N noisy cycle r — 0.1, /3 — 0.5 
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Figure 6. Stochastic models of futile cycles. Empirical probability density function for So at 

t = 10 after 1000 simulations for the futile cycle models with the parameter configurations considered in 

Figure 

fl -1 0\ ai(t) = £(t)[K/(K + P 2 )] 2 <J2 = SrRi 

1 -1 a 3 (t) = {(t)[K/(K + Pi)] 2 a 4 = 5 R R 2 

00001 -1 00 a 5 = apRi a 6 = 5 P Pi 
\0 1 -1/ a 7 = apR 2 a s = 5 P P 2 

Table 4. The bistable model of gene expression in |73) : the stoichiometry matrix (rows in order R\ , 
R2, -Pi) Pz) and the propensity functions. 

bottom panels the cycle with bounded noise and autocorrelation t = 0.1 in (left) and r = 1 in (right). 
In both cases in the top panel the noise intensity is (3 = 0.5 (top) and (3 = 1 (bottom). The initial 
configuration is always (So, E, ESq, Si, F, FS\,N) = (0, 20, 0, 2000, 50, 0, 10) and the kinetic parameters 
are ki = 40, fe_i = k 2 = 10000, k 3 = 200, /c_ 3 = 100, fc 4 = 5000 for the noise-free and the bounded 
noise case, and k$ = ke = 10, fc_g = 5 and fc_6 = 0.2 [22]. We also show in Figure |6] the empirical 
probability density function for the concentration of So, i-e. "P(X(i) = x) given the considered initial 
configuration, at t € {2, 5, 7, 10} after 1000 simulations for the futile cycle models with the parameter 
configurations considered in Figure [SJ The analysis of such distributions outline that for the noise-free 
system the distributions are clearly unimodal, whereas for noisy futile cycle, in both cases, they are bi- 
modal. Moreover, it is important to notice that the smallest peak of the distribution, i.e. the rightmost, 
has a bigger variance when N is considered, rather than when a bounded noise is considered. 

Bistable kinetics of gene expression 

Let us consider a model by Zhdanov [T3J[73] where two genes G\ and G2 , two RNAs i?i and R2 and two 
proteins Pi and Pi arc considered. In such a model synthesis and degradation correspond to 

Gi ^ Gi + Ri Ri^ R1 + P1 R!^* Pi ^ * 
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Figure 7. Periodically perturbed toggle switch. In the top panels a single run for Zhdanov model 
(|33f with a = 0.5 (left) and a = 1 (right). In bottom panels averages of 1000 simulations. In all cases 
ap = 100 mm -1 , otp = 10 mm -1 , 6p = 6p = lmin , K = 100 and r = 100 mi??" 1 and the initial 
configuration is (Pi, Pi, P2, P2) = (10, 0, 0, 0). The populations and the noise are plotted for the single 
run. 



Such a reaction scheme is a genetic toggle switch if the formation of Pi and P2 is suppressed by P2 and 
Pi , respectively [71 fTn[7DH7^] . Zhdanov further simplifies the schema by considering kinetically equivalent 
genes, and by assuming that the mRNA synthesis occurs only if 2 regulatory sites of either Pi or P2 are 
free. The deterministic model of the simplified switch when synthesis is perturbed is 



Ri = f(t) , 
P/ = a P Pi -SpPi 



K 



SpRi 



^ y> \K + Pi 
P 2 ' = apP 2 - SpP 2 



- S R R 2 



(33) 



where the perturbation is 



1 



a sin 



2irt 



Here ap, Sp, ap and Sp are the rate constants of the reactions involved, term \Kj(K + Pi)] 2 is the 
probability that 2 regulatory sites are free and K is the association constant for protein P. Notice that 
here perturbations are given in terms of a time-dependent kinetic function for synthesis, rather than a 
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Figure 8. Periodically perturbed toggle switch. Empirical probability density function at 
t G {900, 950, 1000}, after 1000 simulations for Zhdanov model with the parameter configurations 
considered in Figure [Jj In left a = 0.5 and in right a = 1. 



stochastic differential equation. Before introducing a realistic noise in spite of a perturbation we perform 
some analysis of this model. As in [73] we re-setted model (J33J in a stochastic framework by defining 
the reactions described in Table H] Notice that in there two reactions have a time-dependent propensity 
function, i.e. ai(t) and 03(4) modeling synthesis. 

In the top panels of Figure UJ we show single runs for Zhdanov model where simulations are performed 
with the exact SSA with time-dependent propensity functional. We considered an initial configuration with 
only 10 RNAs Ri. As in [73] we set ap = 100 min -1 , ap = 10 mire -1 , 5r = Sp = lrreire -1 , K = 100 
and t = 100 mire - ; notice that this parameters are realistic since, for instance, protein and mRNA 
degradation usually occur on the minute time-scale [73] ■ We considered two possible noise intensities, i.e. 
a = 0.5 in left and a = 1 in right and, as expected, when a increases the number of switches increases. 
To investigate more in-depth this model we performed 1000 simulations for both the configurations. In 



2 In 1731 an exact SSA 1281 is used to simulated the model under the assumption that variations in the propensity functions 
are slow between two stochastic jumps. This is true for r = 100 as in [73] . but not true in general for small values of r. 
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Figure 9. Periodically perturbed toggle switch. Empirical probability density function for R2 
plotted against time, i.e. the probability of being in any reachable state x for 900 < t < 1000. Lighter 
gradient denotes higher probability values. We used data collected with 1000 simulations of model (|33| 
where r = 100 and two perturbation intensities are used, i.e. a S {0.5, 1} as reported in the top 
captions. In the ir-axis the species concentration is represented, in the y-axis minutes are given. 



the bottom panels of Figure [7] the averages of the simulations are shown. The average of our simulations 
evidences a major expression of protein P 2 against Pi, for both values of a, with dumped oscillations for 
a = 0.5 and almost persistent oscillations for a = 1. 

In Figure |8] we plot the empirical probability density function of the species concentrations, i.e. 
V[X.(t) = x] given the considered initial configuration, at t G {900,950, 1000} min as obtained by 1000 
simulations. Interestingly, these bi-modal probability distributions immediately evidence the presence of 
stochastic bifurcations in the more expressed populations R2 and P 2 . In addition, the distributions for 
the protein seem to oscillate with period around 100, i.e. for a = 1 they arc unimodal at t <E {900, 1000} 
and bi-modal at t = 950. 

For the sake of confirming this hypothesis in Figure [9] the probability density function of P2 is plotted 
against time, i.e. the probability of being in state x at time t, for any reachable state x and time 
900 < t < 1000. In there we plot a heatmap with time on the y-axis and protein concentration on the 
x-axis; in the figure the lighter gradient denotes higher probability values. Clearly, this figure shows the 
oscillatory behavior of the probability distributions for both value of a and, more important, explains 
the uni-modality of the distribution at t = 900 and t = 1000 with a = 1, i.e. the higher variance of the 
rightmost peak at a = 1 makes the two modes collapse. Finally, we omit to show but, as one should 
expect, the oscillations of the probability distribution, which are caused by the presence of a sinusoidal 
perturbation in the parameters, are present and periodic over all the time window < t < 1000. 



Bounded noises. Wc investigated the effect of a Sine- Wiener noise affecting protein synthesis rather 
than a perturbation, i.e. a new £(t) is considered 



£sw(*) = a R 



1 



a sm 




with W a Wiener process. Here simulations are performed by using the SSAn where the reactions in 
Tableware left unchanged, and the propensity functions ai(t) and 03 (t) are modified to 



a 1 {t) = ^ vl {t)[K/(K + P 2 )f 



a 3 (t)=£ sw (t)[K/(K + P 1 )} 
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Figure 10. Stochastically perturbed toggle switch. In top panels, single run for Zhdanov model 
with Sine- Wiener bounded noise: in left a = 0.5, in right a = 1. In bottom panels the averages of 1000 
simulations. In all cases ap = 100 min , ap — 10 min , 5r = 5p = lmin~ l , K = 100 and 
t = 100 mm -1 and the initial configuration is (i?i, P2, R2, P2) — (10,0,0,0). The populations and the 
noise are plotted for the single runs. 



For the sake of comparing the simulations with those in FiguresJTHHl we used the same initial condition 
and the same values for a p., ap, dp, Sp and K. To compare the effect of a realistic noise against the 
original perturbation we simulated the system with the same values i.e. the noise intensity a = 0.5 in left 
and a = 1 in right of the top panels in FigureUHl and in both cases r = 100. As expected, in this case the 
trajectories are more scattered than those in Figure [3 and the switches are still present. However, for 
maximum noise intensity a = 1 time-slots emerge where the stochastic systems predicts a more complex 
outcome of the interaction. In fact, for t G [0,200] U [800,900] neither protein Pi nor P2 seem to be 
as expressed as in the other portions of the simulation, thus suggesting the presence of noise-induced 
equilibria absent when periodic perturbations are present. 

To investigate more in-depth this hypothesis we again performed 1000 simulations for both the config- 
urations, the averages of which are shown in the bottom panels of Figure fTOl In this case, the simulation 
times, which again depend on the noise correlation, span from 3-^5 min to 30 -j- 40 min, thus making 
the choice of good parameters crucial. Differently from the case in which a sinusoidal perturbation is 
considered, i.e. Figure [3 in this case the averages are not oscillatory, but instead show a stable conver- 
gence Also, the final outcome seems again to predict the expression of P2 inhibiting Pi. To understand 
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Figure 11. Stochastically perturbed toggle switch. Empirical probability density function at 
t = 1000, after 1000 simulations for Zhdanov model with Sine- Wiener noise. Parameters are as in 
Figure [TU] and two perturbation intensities are used, i.e. a G {0.5, 1} as reported in the top captions. 



better this point we plotted in Figure [TT] the probability density of reachable states at t — 1000 min, 
i.e. V[X.(t) = x] given the considered initial configuration, and in Figure [12] we plotted that distribution 
against time for R 2 . It is worth noting that we also ranged t over [900, 1000] but since V[X(t) = x] did not 
change we omitted to plot it here. Again, Figure [T^] is a heatmap where on the y-axis time in minutes is 
given, on the x-axis the possible concentration for R 2 and the lighter gradient denotes higher probability 
values. Notice that in this case Figure [T^] represents an empirical evaluation of the solution the DCKE for 
this system, i.e. equation ([TT]) . Both graphics are obtained by 1000 simulations with a = 0.5 (left panels) 
and a = 1 (right panels). These figures show that a low- intensity noise makes the probability distribution 
become three- modal, i.e. notice the two rightmost peaks in Figure [TT] and the white/light-blue gradients 
in Figure [12l Differently, when the noise intensity is higher, the two rightmost peaks almost merge, thus 
forming a bi-modal distribution where the smaller peak almost spreads uniformly on the state space for 
the variables. Notice that, in this case, the amplitude of such a peak is higher than for a = 0.5, i.e. notice 
the intensity of the blue gradient in Figure [T^] For a = 0.5 it is possible to notice two red gradients: one 
approximatively for x — > 200 and one for x G (10, 30). The major peaks in the distribution for R 2 are for 
x < 10, for x <G (50, 100) and for x 6 [130, 180]. The probability of each of these peaks is decreasing as 
x increases, thus confirming the intuition of Figure 1111 Similar considerations can be done when a — 1 
where, as shown by Figure ITTT the first dark-red area separating the first two peaks in a = 0.5 is vanished, 
thus forming a bi-modal instead that a three-modal probability distribution. 

Finally, for the sake of considering a wide range of biologically meaningful values for r, which we recall 
it represents a measure of the speed of noise variation, we evaluated the solution of the DCKE for R 2 for 
the same configuration used in Figure [T2l and r G {1, 10,25, 100} mm. We performed 1000 simulations 
of the model for each value of r with a = 0.5, the value showing a more interesting behavior. In Figure 
[13] the probability of the reachable states at t = 1000 min is plotted. If is immediate to notice that the 
height of the first peak increases as r decreases, and more precisely the distribution seems to switch from 
a three-modal one to a bi-modal when r < 25. In each panel of Figure [TU we plot the variation of such 
probability distribution for 900 < t < 1000. By that figure it is possible to observe that by ranging r 
the dark-red gradient increases in size as far as r decreases. This means that the amplitude between the 
peaks of the density strictly depends on the value of r, thus suggesting a strong role for extrinsic noise 
in determining the network functionalities. 
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Figure 12. Stochastically perturbed toggle switch. Empirical probability density function for i?2 
plotted against time, i.e. the DCKE solution for R2 in 900 < t < 1000. Lighter gradient denotes higher 
probability values. We used data collected with 1000 simulations of Zhdanov model with Sine- Wiener 
noise where r = 100 and two perturbation intensities are used, i.e. a G {0.5, 1} as reported in the top 
captions. In the x-axis the species concentration is represented, in the y-axis minutes are given. 
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Figure 13. Stochastically perturbed toggle switch. Empirical probability density function at 
t = 1000 for i?2, after 1000 simulations for Zhdanov model with Sine- Wiener noise. From left to right 
t = 1, t = 10, t = 25 and r = 100. In all cases a = 0.5 and other parameters are as in Figure ITOl 



T = 100 




50 200 



5C 100 150 



Figure 14. Stochastically perturbed toggle switch. Empirical probability density function for R2 
plotted against time, i.e. the DCKE solution for R2 in 900 < t < 1000. Lighter gradient denotes higher 
probability values. We used data collected with 1000 simulations of Zhdanov model with Sine- Wiener 
noise. From left to right t = 1, r = 10, r = 25 and r = 100. In all cases the noise intensity is a = 0.5. 
In the x-axis the species concentration is represented, in the y-axis minutes are given. 
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Discussions 

In this paper we investigated the effects of joint extrinsic and intrinsic randomness in nonlinear genetic 
networks, under the assumption of non-gaussian bounded external perturbations. Our applications have 
shown that the combination of both intrinsic and extrinsic noise-related phenomena may have a con- 
structive functional role also when the extrinsic noise is bounded. This is in line with other researches 
- only focusing on either intrinsic or extrinsic noise - recasting the classical interpretation of noise as a 
disturbance more or less obfuscating the real behavior of a network. 

This work required the combination of two well-known frameworks, often used to separately describe 
biological systems. We combined the theory of stochastic chemically reacting systems developed by 
Gillespie with Langevin systems describing the bounded variations of kinetic parameters. The former shall 
allow considering the inherent stochastic fluctuations of small numbers of interacting entities, often called 
intrinsic noise, and clearly opposed to classical deterministic models based on differential equations. The 
latter permits to consider the influence of bounded extrinsic noises. These noises are modeled as stochastic 
differential equations. For these kind of systems, although an analytical characterization is unlikely to 
be feasible, we were able to derive a differential Chapman-Kolgomorov equation (DCKE) describing the 
probability of the system to occupy each one of a set of states. Then, in order to analyze these models by 
sampling from this equation we defined an extension of the Gillespie's Stochastic Simulation Algorithm 
(SSA) with a state-dependent Langevin system affecting the model jump rates. This algorithm, despite 
being more costly than the classical Gillespie's SSA, allows for the exact simulation of these doubly 
stochastic systems. 

We outlined the role of extrinsic noise for some biological networks of interest. In particular, we 
were able to extend classical results on the validity of the Michaelis-Menten approximation to the pro- 
totypical Enzymc-Substrate-Product enzymatic reaction by drawing a Stochastic Quasi Steady State 
Assumption (SQSSA) for noisy reactions. Along the line of the classical deterministic or stochastic uses 
of the Michaelis-Menten approximation, this should permit to reduce the size of more general enzymatic 
networks even in presence of extrinsic bounded noises. 

Moreover, we showed that in a recurrent pattern of genetic and enzymatic networks, i.e. the futile 
cycle, the presence of extrinsic noises induces the switching from a monomodal probability density (in 
absence of external perturbations) to a multimodal density. 

Similarly, in the case of the toggle switch, which is inherently multistable, the presence of extrinsic 
noise significantly modulates the probability density of the genes concentration. In this important network 
motif we also investigated the role of periodic perturbations against a realistic noise. 

Thus in general the co-presence of both intrinsic stochasticity and bounded extrinsic random pertur- 
bations might suggest the presence of possibly unknown functional roles for noise for these and other 
networks. The described noise-induced phenomena are shown to be strongly related to physical charac- 
teristics of the extrinsic noise such as the noise amplitude and its autocorrelation time. 

A relevant issue that we are going to investigate in the next future is the role of the specific extrinsic 
bounded perturbations. Indeed, in non-genetic systems affected by bounded noises it has been shown 
that the effects of the perturbations depend not only on the above general characteristics of the noise, 
but also on its whole model [40H421 I76 ] . In other words the transitions of a system perturbed by a sine- 
Wiener noise might be quite different from those induced by another bounded perturbation, for example 
the Cai-Lin noise [77] or the Tsallis noise [33], also when their amplitude and autocorrelation times are 
equal. In other words, a single network in two different environments might show two different behaviors 
depending of fine details of the kind of perturbations that are present. This might also suggest that a 
same network might exhibit many different functions depending on its "locations" . 

Finally, concerning these points, we stress that these peculiar properties of bounded extrinsic per- 
turbations make it even more important the investigations, such as those of [37], aimed at inferring by 
deconvolution the external noise from the experimental data, in order to infer which kind of noise affect 
a given network in a well determined environment. 
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